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I. INTRODUCTION 



Calculations of bulk observables for medium to heavy nuclei typically rely on nonrelativis- 
tic (Skyrme) or covariant (QHD) mean-field models. How can one systematically go beyond 
such models? One possibility is to view nuclear mean-field approaches as approximate im- 
plementations of Kohn-Sham density functional theory (DFT) [|1], ^, ^, |^ , which is widely 
used in condensed matter and quantum chemistry applications. To date, the refinement of 
DFT methods has focused almost exclusively on Coulomb systems and DFT has had little 
explicit impact on nuclear structure phenomenology [0. Effective field theory (EFT), how- 
ever, offers a systematic approach for describing low-energy nuclear physics that can provide 
a framework for nuclear DFT. EFT approaches have been making steady progress on two- 
and three-body nuclear systems and certain halo nuclei, and ab initio shell model methods 



should be able to extend the calculations to a wide range of light nuclei 0, ^, ^, |TD|, |TT 
Our ultimate goal is to systematically describe heavier nuclei by applying EFT methods to 
Kohn-Sham DFT. In this paper we take the first steps toward this goal. 

Density functional theory provides a calculational framework that greatly extends the 
range of many-body calculations in finite systems. In Kohn-Sham DFT, the coordinate-space 
ground-state density for A fermions is found by simply summing the squared wavefunctions 
for a set of single-particle orbitals. These wavefunctions are solutions to a Schrodinger 
equation featuring the Kohn-Sham single-particle potential, which is local and energy in- 
dependent. At zero temperature, for nondegenerate ground states (e.g., closed shells), the 
sum for the density is over the orbitals with the lowest A eigenvalues, equally weighted 
(i.e., "occupation numbers" are either one or zero). The Kohn-Sham potential is itself a 
functional of the orbitals, so we have a self-consistent problem that is solvable by iteration. 
When the iterations have converged, the orbitals can be plugged into an energy functional 
to find the ground-state energy. While this procedure is very similar to the self-consistent 
Hartree approximation, it incorporates all correlations if the correct Kohn-Sham functional 
is used. 

The advantages of the Kohn-Sham DFT procedure are clear. One is that the external 
potential, if present, appears in a very simple way in a separate term in the functional, with 
the rest of the functional being universal. Another is that solving the Schrodinger equa- 
tion with local potentials for eigenvalues is relatively simple and fast. The computational 
advantages would be lost, of course, if the construction and evaluation of the Kohn-Sham 
potential is itself too difficult or expensive. The experience for Coulomb systems is that 
local density approximations (LDA) for the potential part of the energy functional and their 
improvement with gradient expansions work very well in many (although not all) systems 
1^, |12|, |13|. That is, the most important nonlocality to treat explicitly is the kinetic energy, 
which is why the Kohn-Sham framework is generally superior to the Thomas- Fermi approxi- 
mation, in which the kinetic energy is given by an LDA (see Sect. IV D| below). The gradient 
expansion approximations greatly simplify the evaluation of the Kohn-Sham potential. We 
anticipate that the convergence of derivative expansions for nuclear systems will be rapid in 
general, but this will need to be confirmed. 

Our strategy is to apply an effective action formalism H^, |T5|, ^ to calculate the Kohn- 
Sham potential and energy functional order-by-order in an EFT expansion, and to use EFT 
power counting to organize and justify a derivative expansion of the functional. In the 
present work, we use a dilute, confined Fermi system with short-range interactions as a lab- 
oratory to explore how EFT can be used to carry out systematic DFT calculations. In this 
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case, the explicit expansion parameter is the local Fermi momentum times the scattering 
length (and other effective range parameters). We assume a gradient expansion parame- 
ter that justifies a local density approximation, but the verification of this assumption is 
postponed to future work. Ultimately we are interested in calculating self-bound systems 
(e.g., nuclei), with spin- and isospin-dependent interactions and long-range forces (e.g., pion 
exchange). These are all significant but well-defined extensions of the model described here. 
In the meantime, the model provides a prototype for more complex systems and also has 
a physical realization in recent and forthcoming experiments on fermionic atoms in optical 
traps fT^j. 

The Kohn-Sham approach to DFT was proposed in Ref. Since then, the literature of 
DFT applications has grown exponentially, primarily in the areas of quantum chemistry and 
electronic structure [§. A general introduction to density functional theory as conventionally 
applied is provided in the books by Dreizler and Gross and Parr and Young 0, while 



Ref. [0 is a practitioners guide to DFT for quantum chemists. The connection of DFT 
to nonrelativistic mean-field approaches to nuclei (e.g., Skyrme models) was pointed out 
in Ref. |jT9| (and no doubt elsewhere) and was explored for covariant nuclear mean-field 



models in Refs. pO|, |21|. However, it has not led, to our knowledge, to new or systematically 
improved mean-field-type functionals for nuclei. 

The use of functional Legendre transformations for DFT with the effective action formal- 



ism was first detailed by Fukuda and collaborators p2|, |23|, who also discuss the inversion 
and auxiliary field methods of constructing the effective action. The connection to Kohn- 
Sham DFT was shown by Valiev and Fernando |]2^, ^ |2^, ^ and later by other authors in 



Refs. [p8| , p9| , [30[] . Recent work by Polonyi and Sailer applies renormalization group methods 
and a cluster expansion to an effective-action formulation of generalized DFT for Coulomb 
systems To our knowledge, however, there is no prior work on merging the Kohn-Sham 



density functional approach and effective field theory. 

The plan of the paper is as follows. In Sect. we review effective field theory for a 
dilute system of fermions. In Sect. [in| , the effective action approach for determining a 
Kohn-Sham energy functional through a Legendre transformation is reviewed. A systematic 
approximation procedure for constructing the energy functional, the inversion method, is 
presented. In Sect. the formalism is applied to a dilute Fermi gas in a harmonic trap 
and results are presented through third order (NNLO) in the dilute EFT expansion using 
an LDA. Section summarizes our results and future plans. 

II. EFT FOR INFINITE DILUTE FERMI SYSTEMS 
A. Background 

Effective field theory (EFT) provides a powerful framework to study low-energy phenom- 
ena in a model-independent way ||^, ^ 32|. The EFT approach is grounded in some very 



general physical principles |^ . If a system is probed or interacts at low energies, resolution 
is also low, and fine details of what happens at short distances or in high-energy interme- 
diate states are not resolved. Therefore, the short- distance structure can be replaced by 
something simpler without distorting the low-energy observables. This is analogous to a 
multipole expansion, in which a complicated, charge or current distribution is replaced for 
long-wavelength probes by a series of point multipoles. EFT uses local Lagrangian field 
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theory as a framework for carrying out this program in a complete and systematic way.^ 
The uncertainty principle implies that high-energy intermediate states are highly virtual 
and only last for a short time, so their effects are not distinguishable from those of local 
operators ||3^. This physics can then be systematically absorbed into the coefficients of 
these operators using renormalization. 

The effective degrees of freedom (dof 's) in an EFT depend on a separation or resolution 
momentum scale A, which sets the radius of convergence of an EFT. Long-range dof's with 
respect to A must be treated explicitly while short-range physics is encoded in the coefficients 
of the local operators. (See Ref. for details of how the pion can be considered either a 
short- or long-distance degree of freedom in the two-nucleon problem, depending on the 
resolution scale.) The hierarchy of scales in the system is exploited to provide expansion 
parameters. For example, if the typical momenta k are small compared to the inverse range 
of the interaction 1/R, we can take A ~ and low-energy observables can be described by 
a controlled expansion in kR. All short-distance effects are systematically absorbed into low- 
energy constants through renormalization. As a result, the EFT approach allows for accurate 
calculations of low-energy processes and properties with well-defined error estimates (based 
on the order of truncation) 0, ||, |^, |T0|. 

The application of EFT methods to many-body problems promises a consistent organi- 
zation of many-body corrections, with reliable error estimates, and insight into the analytic 
structure of observables (see, for example, the identification using renormalization group 
methods of logarithmic contributions to the energy of dilute systems in Refs. |Q and [0). 
The EFT provides a model-independent description of finite-density observables in terms of 
parameters that can be fixed from scattering in the vacuum or from a subset of finite density 
properties. One can also exploit the freedom in an EFT of using different regulators and 
renormalization schemes to find simplifications and clarifications [p^ . 

While EFT has shown early promise in applications to basic many-body problems (e.g., 
Refs. and [^), there are formidable challenges in carrying out many-body calculations, 
particularly for finite, non-uniform systems. For sufficiently small numbers of fermions, the 
many-body Schrodinger equation for finite systems can be solved directly, for example by 
Green's function Monte Carlo methods. However, the computational cost of these methods 
grows as a power of the number of particles (or faster), which prevents their application 
to very many systems of interest in condensed matter and quantum chemistry (Coulomb 
systems) and nuclear physics (medium to heavy nuclei). 

The purpose of the present work is to address these challenges systematically. We seek 
to merge the organizational advantages and insight provided by EFT with the calculational 
power and relative ease of DFT for finite systems. We will build the basic EFT formalism 
of this merger for a dilute gas of identical fermions with short-range interactions. We re- 
view below the results from the calculation of the energy density in the case of a uniform 
system [^]. These results will be the starting point of the DFT calculation of the energy 
density of a dilute Fermi gas confined by an external harmonic potential. They will also 
serve as the basis for the LDA to the energy functional. 



Note that conventional nuclear phenomenology also relies on these principles in using potentials cut 
off at short distances. However, cutoff independence is a goal of EFT that is not usually achieved in 
phenomenological approaches. 
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B. Lagrangian and Energy Density for Uniform System 



We consider a general local Lagrangian for a nonrelativistic fermion field that is invariant 
under Galilean, parity, and time-reversal transformations: 

L = + + — ]v^-^(V^V)2 + ^[(V;^)t(V;V2^)+ h.c. 

+ + , (1) 

where V is the Galilean invariant derivative and h.c. denotes the Hermitian 

conjugate. The terms proportional to C2 and C2 contribute to s-wave and p-wave scattering, 
respectively, while the dots represent terms with more derivatives and/or more fields. The 
Lagrangian Eq. (|l|) represents a particular canonical form, which can be reached via field 
redefinitions. For example, higher-order terms with time derivatives are omitted, as they 



can be eliminated in favor of terms with spatial derivatives |^ . 

To reproduce the results in Ref. ||3^, we can write a conventional generating functional 
with the Lagrangian of Eq. (|1|) and Grassmann sources coupled to ip'^ and ip, respectively 



| 35| . The non-quadratic part of the Lagrangian is removed in favor of functional derivatives 



with respect to the Grassmann sources and the remaining quadratic part is evaluated in 
terms of a non-interacting Green's function times the sources. Perturbative expansions 
for Green's functions (and subsequently S-matrix elements) follow by taking successive 
functional derivatives, and the ground state energy density follows by applying the linked 
cluster theorem (see Ref. for details). In calculating the energy, finite density boundary 
conditions at T = can be incorporated into the non-interacting Green's function using the 
chemical potential fi or by including them by hand with a non-interacting chemical potential. 



The latter approach is simplest at T = and was adopted in Ref. [0. 

The coefficients Cq, C2, and can be obtained from matching the EFT to a more 
fundamental theory or to (at least) three independent pieces of experimental data. We 
follow the regularization and renormalization prescription described in Ref. [Q, namely 
dimensional regularization with minimal subtraction, which is particularly convenient for 
the dilute, natural system. By matching to the effective-range expansion for low-energy 
fermion-fermion scattering, we can express the C2i in terms of the effective-range parameters: 

C„ = ^^, C. = C„^, and C^^'-^. (2) 

where (op) are the s-wave (p-wave) scattering length and is the s-wave effective range, 
respectively. 

The energy density for the uniform dilute fermion system with natural scattering length 
33| is calculated as a perturbative expansion in kp/A, where k-p is the Fermi momentum 



and A is the resolution scale (e.g., A ^ 1/R for hard spheres). The non- interacting energy 
density at zero temperature for A particles with spin-degeneracy g in volume V can be 
written as ^ 



where the density p is 
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FIG. 1: Hugenholtz diagrams for a dilute Fermi gas through order kp in the energy density. 



The order-by-order corrections to Eq. (^) due to interactions can be represented by the 
Hugenholtz diagrams given in Fig. |I|. The Feynman rules for calculating these graphs along 
with details of the renormalization needed to render divergent graphs finite are given in 
Ref. ^ 
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Here we will only need the final results and we simply quote them up to next-to- 
next-to-leading order (NNLO). 

The LO diagram of order kp [Fig. 0(a)], which represents the Hartree-Fock result, con- 



tributes 



Si = p{g-1] 



kp 



kptts 



(5) 
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where Eqs. (0) and (^) have been applied. At order kp (NLO) there are two diagrams. The 
three-loop diagram [Fig. 0(c)] is an example of an "anomalous diagram," which vanishes 
identically for a uniform system in the zero-temperature formalism but is nonzero when cal- 
culated in the zero-temperature limit. In the latter case, its contribution is precisely canceled 
by the shift between the noninteracting and interacting chemical potentials, as dictated by 



the Kohn-Luttinger-Ward theorem |^5|, |3^. We will find an analogous cancellation in 
the diagrammatic expansion of the Kohn-Sham functional. The other diagram at this order 
[the "beach ball" diagram of Fig. |T](b)], makes the contribution 



S2 = p{g-1] 



kp 



(kptts 



'11 -21n2l 



(6) 



2M 35772 

Finally, we have the graphs of order kp (NNLO). The first three are anomalous diagrams. 
In addition to graphs containing the Cq vertex [Figs. |I](g) and (h)], there also graphs that 
contain C2 and Cg. Altogether these graphs give the correction 



^3 



P 



kp 
2M 



(9-^)Y^ikpa 



' kpTs + {9 + ^) ^ {kpap) 

OTT 



+ {g - 1){(0.07550 ± 0.00003) + (^ - 3) (0.05741 ± 0.00002)} (kpa. 



(7) 



where the integrals for Figs. |I](g) and (h) have been evaluated numerically. The expansion 
can be continued systematically, but this is as far as we will need here. 
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III. EFFECTIVE ACTIONS AND THE INVERSION METHOD 



The energy density given in the last section illustrates that the ground state energy of 
a uniform, interacting many-body system may be written as a function of the constant 
density. The DFT theorems of Hohenberg and Kohn (HK) formally prove the existence 
of a generalization to finite, nonuniform systems: the ground state energy can be obtained 
from a functional of the local density alone. Specifically, for a system with an external 
potential f (x) coupled to the density p(x), there exists an energy functional E[p] that can 
be decomposed as 

E[p(x)] = Fhx[p(x)] + J rf=^xt;(x)p(x) , (8) 

where the functional Fhk[p], which is known as the HK free energy, is universal (i.e., in- 
dependent of the potential v). A variational principle ensures that the functional E[p] is a 
minimum equal to the ground state energy when evaluated at the exact ground state den- 
sity. The practical problem of density functional theory is to find accurate and tractable 
approximations to Fhk[p] 0; 0; 3- 

In applications to Coulomb systems, Fhk[p] is typically decomposed into a noninteracting 
kinetic energy Fni[p], the Hartree term Eh[p], and everything else, which is defined as the 
exchange-correlation energy E^dp] 0: 

Fhk [p] = [p] + Eh [p] + E^, [p] . (9) 

In the Thomas-Fermi approximation, the noninteracting functional is evaluated within an 
LDA. This approximation is accurate at high densities when the characteristic scale over 
which the fermion density changes is small compared to the Fermi wavelength. The approx- 
imation is inadequate for most ordinary Coulomb systems and improves only when some 
account is taken of the nonlocality in Fni- Kohn and Sham [|l| introduced a method involving 
auxiliary orbits by which Fni could be treated exactly, which has become the basis for most 
practical DFT applications 0, |^ . 

The utility of DFT then rests on finding an explicit expression, exact or approximate, for 
Exc- For Coulomb systems, the conventional procedure is to approximate E^c in the LDA 
based on a fit to Monte Carlo results for the energy of a uniform electron gas as a function 
of the density, and then to include semi-phenomenological gradient corrections [^]. For 
applications to other systems, such as nuclei or trapped atoms, we adopt an approximation 



scheme based on effective action methods. Here we follow the work of Fukuda et al. [0, 
and its further development by Valiev and Fernando [^. We extend this work by merging 
it with an effective field theory expansion. For consistency with the EFT treatment of dilute 
systems in Ref. [^, we work at zero temperature in Minkowski space. 

We begin with the system described by the Lagrangian in Eq. (p, to which we add a 
term for an external potential t'(x) coupled to the density operator, v{^iIj'^iI)? Here we take 



^ In general, the potential will be coupled to more operators than simply 'ip'^ip. Field redefinitions, which 
leave observables unchanged, could be used to eliminate couplings to these operators, which would be 
higher order in the EFT expansion [^9[ . However, such field transformations would also induce energy- 
dependent terms in the Lagrangian. Since we have already used field transformations to achieve a canon- 
ical, energy-independent Lagrangian, we will assume that omitted higher-order couplings, which should 
be suppressed numerically according to power counting, can be treated perturbatively. 
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the external potential to be an isotropic harmonic confining potential 

1 



t;(x) = -mcu^ |x|% (10) 



as might be appropriate for some atomic traps but the discussion holds for a general 



external potential. We also introduce a c- number source, J{x), coupled to the composite 
density operator and write down the generating functional using the path integral formula- 
tion, 

Z[J] = e'^t-'l = j DtfjD^^ e'^'''^ [£ + JWt(,)^(,)] _ ^^^^ 

For simplicity, normalization factors are considered to be implicit in the functional integra- 
tion measure. (See Refs. |2^, for a more careful treatment of the path integrals.) Using 



the definition in Eq. (IT^), we see that the density (in the presence of J) is 

p{x) ^ = . (12) 

dJ[x) 

While it is possible to absorb f (x) into the definition of J{x), we find it more convenient to 
recover our original system in the limit J — ^ 0. 

The effective action is defined through the functional Legendre transformation 

r[p] = W{J] - J (fx J{x)p{x) . (13) 

This transformation ensures that F has no explicit dependence on J. The existence of such a 
functional follows from the concavity of W[J] which guarantees that Eq. (|12D can be inverted 
to give J = J[p\. The functional dependence between J and p can be used to write Eq. (|13D 
entirely in terms of p. The proof that W[J] is strictly concave is given in Ref. [2^ . 

Since we are interested here in describing finite many-body ground states at T = 0, it is 
most convenient to work with functions of the particle number A rather than the chemical 
potential p. This can be achieved via a conventional Legendre transformation on W or F, 
but is more simply carried out implicitly, by choosing appropriate finite-density boundary 
conditions that enforce a given A by hand (the actual procedure is detailed below). In 
the following, we will assume this has been done. Thus, F and W are functions of A and 
variations over p{x) conserve A. In addition, we restrict the discussion to time independent 
sources and densities. In this case the effective action acquires a factor that corresponds to 
the time interval over which the source is acting, which we indicate schematically as 

/oo 
dt , (14) 
-oo 

as in Ref. |^. To avoid overly cluttered notation, we will divide out this ubiquitous time 
factor everywhere it appears and write 



F[p]=F[p]x 



oo 

dt 

-oo 



-E[p] . (15) 



and similarly with W[J] and the expansions below. (We will continue to use F rather than 
E in this section.) 



8 



In conventional treatments (e.g., see Refs. [|T^, 0), an effective action is derived from 



a Legendre transformation with respect to a source coupled to one of the fields in the 
Lagrangian, rather than to a composite operator as in the present case. However, the usual 
advantages of working with an effective action are also present here. The effective action 
has extrema at the possible quantum ground states of the system, and when evaluated at 
the minimum is proportional (at zero temperature) to the ground state energy |2^, 531, ET 



In particular, Eq. (|I^) defines an energy functional E[p\ equal to the ground-state energy 
when evaluated with the ground-state density p. 

The extremization condition is shown as follows. Combining Eq. ( [T^ ) and Eq. (|TBp we 
find _ 

or ^ 



The invertibility of ([T^) implies 
so we must have 



^ , (18) 



5J(x 

5f[p\ 



-J(x) . (19) 



5p(x) 

The above equation tell us that when J(x) = the effective action is extremized, which is 
a statement of the second HK theorem [^]. The strict concavity of W[J] implies the strict 
concavity of r[p] and so the extremum is a maximum. Since </(x) = corresponds to the 
original system we see that the energy functional is minimized when evaluated at the exact 
expectation value of the density. 

We observe that the separation of a w(x)-dependent part of the DFT energy functional 



(or r) from a universal part follows directly in the effective action formalism ||22[. From the 
definitions of Z and W , it follows that 

W,=q[J] = W[J + v] (20) 

for any J(x). If we designate Jp(x) the inversion of 5W /5J = p and Jp(x) the inversion of 
SWy=o/6J = p for the same density, then 

5J(x) 5J(x) 5J(x) ' ^ 

which implies 

Jp = J'i + v. (22) 
Upon substituting Eq. (p2D into Eq. (0), the effective action becomes 

f[p] = w^[jp]-y"rf^xj,(x)p(x) 

= W.=o[Jp] - I rf'x (j°(x) + v{^)) p(x) 
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iy.=o[Jj] - J rf'x J^x)] - J d^^v{^) p(x) 
r,=o[p] - / rf=^xt;(x)p(x) , (23) 



which is the promised resuh [Eq. (P)] with an overall minus sign. 

Now consider a system that can be characterized by an effective field theory power- 
counting parameter, which we label A. This parameter may be dimensionful and can appear 
to all orders. We will not exhibit it explicitly here but indicate by subscripts the order in A 
of a given function or functional. In previous discussions of the inversion method p3| , 
the parameter A was always a coupling constant (e.g., for the Coulomb interaction). In 
contrast, we associate A with an appropriate EFT expansion parameter. For example, A 
could be 1/A in the dilute expansion (which we use here) or in a large N expansion, 
where 'W is the spin degeneracy g |]3^. The effective action functional will depend on A 
but we treat p and A as independent variables: 

f = f[p,A]. (24) 

However, the ground state expectation value Pg(x) will naturally depend on A as determined 
by 

5f[p,X] 



5p(x) 



= . (25) 

P=Pg 



That is, if A is changed, a different pg will be necessary to satisfy this equation. 
The Legendre transformation defining F is 

f[p,X]=W[J,X]-J{l)p{l) , (26) 
where J is a functional of p and A as well, as dictated by 

Here we have introduced the convenient shorthand notation 

A{1)B{1) = jd'^, A(xi)5(xO . (28) 

As A is changed, J must be adjusted so that the same p is obtained when taking this 
derivative. This is how the dependence of J on A arises; clearly this dependence can become 
quite complicated. 

The inversion method now proceeds by expanding each of the quantities that depend on 
A in Eq. in a Taylor series in A: 

^[p. A] = Mp] + Ji[p] + J2M + ■ ■ ■ , (29) 
W[J, A] = Wo[J] + Wi[J] + W2[J] + --- , (30) 

r[p,A] = ro[p]+fi[p]+f2[p] + --- , (31) 

where, as advertised, the power of A associated with each function or functional is indicated 
by the subscript. We can substitute the expansion for J into the expansion for W and do 
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a functional Taylor expansion of W[J] about Jo; this makes the A dependence manifest. 
Equating equal powers of A gives (/ = 0, 1, 2, . . .) 



k=l 

I T fciH \-km<l 



+ V— T V jr J TTT^^fci 1 ■■•^fe^ m • 32 

m! dJn(l)---oJo(in) 

m=2 fci,--,fe„>l ^ uv ; 



Since p is independent of A each Jfc[p] follows from each F^: 

..(D^-gM. ,33) 

We reiterate that all of the J;'s as defined here are functionals of p. 
Starting with the zeroth order expression, 

fo[p] = Wo[Jo]-Ml)p{l) , (34) 
we take its functional derivative with respect to p: 

Rearranging, 

■sm-M ^(i,)\WO^o. (36) 



which implies 



since the strict concavity of Tq[p] prohibits 6Jo{l')/6p{l) from having zero eigenvalues. 

Equation ( P7| ) says that Jo(x) is the source (or potential, see below) that generates the 
expectation value p from the noninteracting system (that is, the system defined by A = 0, 
which includes the external potential and Jo(x) but no interactions). Jo(x) is not an arbitrary 
function, but the particular one that has this property. The existence of a Jo{x.) with this 
property is the cornerstone of the Kohn-Sham formalism. 

Equation (|37D also implies that the second term in Eq. ( ^21) cancels with the k = I term 
of the first sum for all I > 0, and thus Ti simplifies to 



m = W,[Jo]-5,oMl)p{l) + Y.^-^^^Ml) 



k=l 

I ^ fei4---+fcm</ 



+ 1^^, L SMl)---5Mm) J^^^^)---J^^^^)- (38) 

m=2 fei,---,fc,„>l ' UV / 

These equations allow us to build the F^'s recursively. Note that the functionals have 
the same diagrammatic expansion as in Fig. ^ but the fermion lines are evaluated with 
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Kohn-Sham (KS) propagators (see below). For a given /, we only need W^s with k less 
than or equal to I and J^'s with k smaller than I (which means lower-order F^-'s). We will 
illustrate the procedure by constructing the first few orders. 

Since the lowest-order term in W[J] corresponds to the system without interactions be- 
tween the fermions, we can write VFo[<^] explicitly by introducing normalized single-particle 
orbitals that satisfy the equation 

+ v{x) - Jo(x) (fi{x) = £i<^i(x) . (39) 



2M 



The index i represents all quantum numbers except for the spin (we consider only spin- 
independent interactions here). Wo[<^o] is then (minus) the sum of the single-particle eigen- 
values up to the Fermi energy ep (which is equal to the chemical potential) 

WolM = ~gY,e,. (40) 

ei<eF 

[Equation ( ^OD can be derived by evaluating Wo[-^o] oc Trln(G°g)~^ using the Kohn-Sham 
Green's function defined below in Eq. (|48|).] In practice, e-p is determined by simply 
counting orbitals until the A lowest are filled (accounting for the spin degeneracy g). Note 
that the e^'s are functionals of Jo through Eq. (|39|); using the normalization of the yji's, we 
find 



j d'x^li^) + v{x) - Jo(x)^ (/..(x) = -^*(y)y..(y) . (41) 



^^o(y) SMy) 

Equations (^7|) and (^Oj) show that the density may be written as 



occ. r OCC. 

V '^^o(y) i 

where the sum is over occupied ("occ") states. Equation (^^ corresponds to the famous 
result of Kohn and Sham, which gives the exact ground state density in terms of the orbitals 
of a non- interacting system |Q. 

With the above results, the lowest order effective action is 

occ. „ 

ro[p] = -gY,^^~ / M^) P(x) . (43) 
We can also use Eq. (|39D to eliminate Ei from Eq. (|^) so that it reads 

ro[p] = -r,[p]-y"rf^xt;(x)p(x), (44) 



where 



occ. „ . 2 \ 

Ts[p]=gJ2j d'x ^*(x) ^-^j ^.(x) (45) 



is the total kinetic energy of the KS non- interacting system. If Fo[p] were given as an explicit 
functional of p, then Jo[p] could be determined by taking a functional derivative according 
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to Eq. (|35|). However, taking the functional derivative of the expression in Eq. ( ^3]) merely 



reproduces the result of Eq. (|37|). Instead, we follow Ref. [2^] and determine Jo from the 
interacting effective action, which we now construct. 
From Eq. (^) we can find ri[p] since 

fM=W^[Jo[p]] . (46) 

For the dilute Fermi system, this is easily calculated. We first introduce the Green's function 
of the KS non- interacting system, x'), which satisfies 

+ ^ - ^(^) + ^o(^)) ^°s(xt, xY) = 6%^ - x')<5(t - t') (47) 



with finite density boundary conditions |^. The KS Green's function has the usual spectral 
decomposition in terms of the orbitals of Eq. (|39|) : 

^Gl{^t, xr) = J2 ^'Kx) ^*(x') e— [9{t - t') d{e, - e^) - d{t' - t) 9{e^ - e,)] . (48) 

i 

The Feynman rules in position space follow conventionally from Eq. (|l]) 0, ^ and we have 
from Fig. |^(a) [with fermion lines representing iG^^ 

Wi[M = ^gig-l)CoJ rf^x <(x, a;+) Gl{x, a;+) , (49) 



where the right side is independent of Xq by Eq. (^H]). The Green's function with equal 
arguments can be directly expressed in terms of the density, 

pi^) = -tgGlix,x^) . (50) 

Using this result and Eq. (^6D, we have 

fiM = -i^^Co I rf^x|p(x)p. (51) 

Since the dependence on p(x) is explicit in Eq. (|5l|) , we can directly take the functional 
derivative with respect to p to obtain 

J,(x) = ^^^^ p(x) . (52) 

Direct functional derivatives with respect to p will not be possible at higher order. How- 
ever, we can find functional derivatives with respect to Jo- An alternative path to Ji[p] from 
fi[p] is 

5p(x) J oJo{y) dp(x) J oJoiy) 
which defines the inverse "density-density" correlator 



^ 5p(x) U^o(x)5Jo(y) 



(54) 
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r,„.= oo+{n)+cxC>o+oOOo 



oo 




FIG. 2: Hugenholtz diagrams for the LO and NLO contributions to the Kohn-Sham interaction 



effective action Tint- The cancellation of the last two diagrams on the first line is given by Eq. (61). 



We can find an expression for D ^ hy taking the functional derivative with respect to Jo in 
Eq. (^91). Remembering that is a functional of Jq it is straightforward to show that 



5JnU) 



Gl^{xi,x)Gl^{x,X2) dxo . 



(55) 



Using Eq. (^5]), the derivative of Wi becomes 

6Wi[J] _ 
5Jo(x) " 



g{g-l) Co Jd'ydxo Gl{y,x)Gl{x,y+)Gl{y,y- 
i{9 -'^) Co d^ydxo Gl,{y , x) Gl,{x , y^) p{y) . 



(56) 



Comparing Eq. (|56|) to Eq. (^) we find 



^-^(x,y) 



9 



dyo Gl^{x,y)Gl^{y,x) 



(57) 



This correlator will appear in all higher-order contributions. 
Having determined Ji[p], we can find T2[p] from 



r^lp] = W2[Jo]+ /rf3x^^^Ji(x) + ± Id^^d'y 



5Jo(x) 5Jo(y 



■^i(x) Ji(y) 



(58) 



= W2[Jo] + - d'^d'y ' D i x,y — f-i. 

2 J dJo{x) oJo{y) 

W2[Jo] is calculated from the graphs Figs. 0(b) and (c): 

W2[Jo] =tg{9-l)^ Jd'x d% Gl{x,y)Gl{x,y)Gl{y,x)Gl{y,x) 

- ig{g -l?^ I d'x d'y G°,(x, x+)G°,(x, y)Gl{y, x)Gl{y, y^) . (59) 



By using Eqs. ([56| ) and (|57D, we can show explicitly that the second term in the expression 
for F2 exactly cancels the second term in W2- First note that 
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FIG. 3: General cancellation of the inverse correlator D for zero-range interactions with no 
derivatives. 



n -1 



X jd'wj dyo Gliw, y)Gl{y, w^)Gl{w, w+) 
-tCo{g-l)GUx,x+). (60) 



We now see that 



2 J 5Jo(x) ^ '"^ 5Jo(y) 

= wig - 1)' T / ^'"^ ' (6^) 

which is the negative of the second term in Eq. (0). This cancellation is shown in Fig. 
where is represented by a double line. Thus, Eq. (|58|) reduces to 

f = ig{g - 1) ^ / ^'2/ Glix, y)Gl{x, y)Gl{y, x)Gl{y, x) . (62) 

The second term in W2 corresponds to the "anomalous" graph (c) in Fig. |l]. This graph 
vanishes identically in the uniform system at zero temperature. However, it does not vanish 
at finite temperature or in a finite system (at any temperature). Rather, the cancellation ex- 
hibited above is analogous to the Kohn-Luttinger-Ward theorem mentioned earlier. Similar 
cancellations, of the type illustrated diagrammatically in Fig. ^, completely eliminate con- 
tributions of to the effective action up to N^LO in the EFT expansion for short-range 
forces. This complete cancellation does not occur for long-range forces, or if the zero-range 
delta functions at the Co vertices are regulated by a cutoff rather than by dimensional 
regularization, as used here. 

We will illustrate the construction of graphically to indicate how the cancellations 
occur at NNLO. We start with the explicit functional expression: 

Figure ^ illustrates that the two terms with J2 factors cancel with each other. Every func- 
tional differentiation with respect to Jo, which inserts a density operator, is represented by 
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FIG. 4: Graphical representation of the Kohn-Sham potential Ji(x) from Eqs. (52) and (|5 



FIG. 5: Cancellation of contributions to involving J2 [see Eq. (p3D]. The graphical representation 
of Ji (with an important minus sign) comes from Eqs. (^) and (^0|). 



a small dot (the large dot is a Cq vertex). VTi corresponds to Fig. |T](a), so the functional 
derivative in Fig. |^(a) yields the diagram shown; an explicit diagram for J2 is not needed. 
Similarly, the second term is represented in Fig. |^(b). As shown, substituting Ji (see Fig. ^ 
yields minus the first contribution, so the sum is zero and J2 does not appear in Fa. 
Thus, F3 reduces to 

W3 is given by the sum of Hugenholtz graphs in Fig. |l|(d) through (j). But the sum of the 
other terms in Eq. (0) exactly cancel the diagrams in Fig. |l](d), (e), and (f), as shown 
in Fig. 1^. The factors in front of the Feynman diagrams indicate additional multiplicative 
factors beyond those prescribed by the Feynman rules. These factors conspire to precisely 
subtract the anomalous diagrams in W3, leaving only Fig. |l|(g) through (j). 

All higher orders in r[p. A] are determined in a similar manner. Direct calculation be- 
comes cumbersome, but one can formulate Feynman rules for F, which dictate how to make 
appropriate insertions of D^^; these are given in Refs. P3| , It is important to 



note that the Kohn-Sham potential Jo completely determines each order in the expansion 
off. 

The cancellations exhibited here at low orders are expected from a comparison of the 
DFT/EFT expansion to perturbation theory. In perturbation theory for the confined sys- 
tem, the energy is calculated by evaluating the diagrams in Fig. |l] using the noninteracting 
propagator in the presence of the external potential for the fermion lines. (This propagator 
would take the form of Eq. (^8[) with harmonic oscillator wave functions for the orbitals.) 
Using the self-consistent Gks instead sums an infinite class of higher-order diagrams at each 
order. The cancellation of anomalous diagrams from F corresponds to the removal of contri- 
butions already included through self-consistency (e.g., tadpoles). Although the nonpertur- 
bative contributions for the dilute, natural system are not required by power counting, the 
self-consistent calculation of the energy and density together is actually easier in practice 
than the purely perturbative calculation. 
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Jl 



+ ^lOOO + oOo 
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OOOO + 3x 
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FIG. 6: Cancellation of contributions to with anomalous diagrams in Fig. |^. The graphical 
representation of Ji (with an important minus sign) comes from Eqs. ( ^2[) and dsT 



The appearance of the inverse density- density correlator D^^ can be understood by com- 
parison to the effective actions for local fields and non-local composite fields. In the former 
case, the Legendre transformation removes one-particle intermediate states (leaving only 
one-particle-irreducible diagrams), while in the latter case, the Legendre transformation re- 
moves two-particle intermediate states (leaving two-particle-irreducible diagrams). Thus we 
infer that the role of is to remove intermediate states created by ip'^ip. The difference 
here is that we cannot write a closed-form expression for the effective action, which is pos- 
sible in the other cases. However, we have seen that for short-range interactions, the extra 
diagrams at low orders cancel against anomalous diagrams, which is a great simplification. 

To find an expression for Jq, we apply the variational principle satisfied by T[p, X] to its 
expansion: 



5T[p,X] 



5p(l) 



P=Pgs 



5p(l) 



-^o(l)L.+ 



P=Ps 



5pil) 



(65) 



p=pg 



where the interaction effective action is 



(66) 



i=l 
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FIG. 7: Contributions to the Kohn-Sham potential Jo through NLO [see Eq. (| 



or 



^o(i: 



p=pg 



5p(l) 



5MV 



P=Pgs 



(67) 



P=Pe 



We stress that these relations hold only when we are solving for the Kohn-Sham potential 
corresponding to the ground-state density Pgs(x). 

The second equality in Eq. (^), which is shown diagrammatically through NLO in Fig. 
is the key to the general Kohn-Sham self-consistent procedure: 

1. Choose an approximation for rint[p] by truncating the expansion in A at some order. 

2. Make a reasonable initial guess for the Kohn-Sham potential Jo(x). 

3. Calculate TintH starting from Jq. 

4. Use Eq. (|67|) to determine a new Kohn-Sham potential Jo(x). 



5. Repeat the last two steps until self-consistency is reached (i.e., until some measure of 
the change in Jq is less than a given tolerance). 

In the next section, we will apply an LDA to Tint, which means that it will given explicitly as 
a functional of p. In this case, the second equality in Eq. (^) is superfluous. One can also 
avoid explicitly evaluating Eq. (0) by adjusting Jo(x) using a steepest-descent approach 

To this point we have neglected to mention that expressions such as r2 in Eq. 



521) are 

divergent. For a uniform system, Jq is constant. In this case, W2 and subsequently r2 are 
renormalized by dimensional regularization and minimal subtraction as in Ref. In a 

finite system, with Jq a function of x, the ultraviolet linear divergence in Eq. (B^) is renor- 



malized by the same counterterm |Q, but it is computationally awkward to renormalize in 
the finite system. By using a derivative expansion, we can perform all renormalizations in 
the uniform system. This will be carried out explicitly in future work. With the LDA trun- 
cation applied below, we can simply use the renormalized expressions for the energy density 
from Ref. to calculate the renormalized Kohn-Sham potential and energy functional. 
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IV. RESULTS FOR DILUTE FERMI SYSTEM IN A TRAP 



In this section, we present representative numerical results for the dilute Fermi system 
defined in Sect. || when confined in a harmonic oscillator trap. Our principal goal here is 
to illustrate the nature of the convergence of the EFT in a finite system both qualitatively 
and quantitatively. We consider effective range parameters corresponding to both attractive 
and repulsive underlying interactions (but we do not allow for pairing). To emphasize the 
difference between Kohn-Sham (KS) and Thomas-Fermi (TF) approaches and since we are 
ultimately interested in nuclear systems, we consider relatively small numbers of trapped 
fermions. Current experiments with trapped atoms use 10^-10^ atoms |T^; as noted below, 
for such large systems the differences between KS and TF get washed out. 



A. The Local Density Approximation (LDA) 

To carry out the Kohn-Sham self-consistent procedure, we need to evaluate the expres- 
sions in the expansion of T[p], so that we can use Eq. ( |67D to find Jo(x). In applications to 
Coulomb systems, the Kohn-Sham energy functional is conventionally written as 

E[p] = T,[p] + j rf=^xt;(x)p(x) + E^[p] + E,,[p] , (68) 

where Ts is given in Eq. (^), -Eh is the Hartree energy and Exc[p] is known as the exchange- 
correlation energy. The Hartree energy is singled out because it can be explicitly written in 
terms of the density. When we have contact interactions only, the Fock term has the same 
dependence on the density as the Hartree term, and so we can also include it explicitly and 
replace E-^ by -Ehf, redefining E^c appropriately. 



This decomposition corresponds to writing the effective action as p4 



r[p] = ro[p] + ri[p] + ^r,,[p] , (69) 



i=2 



up to an overall minus sign. The expression for Fq given in Eq. (H) shows that it has the 
same form as the first two terms in the energy functional in Eq. (|6^). In general, the explicit 
functional dependence of Tg on the density is unknown since it enters implicitly through 
the orbitals. In practice it is easiest to calculate the first two terms by writing them as in 
Eq. (HD: 

f o[p] = -gY,e,- I t/^x Jo(x) p(x) . (70) 

The other terms in F may be written using the Green's function for the KS non-interacting 
system as shown in Sec. |ITT[ The KS Green's functions are explicit functionals of Jo not p, 
however, and therefore almost all of F is not given as an explicit functional of the density; 
the general exception is the Hartree term and here the Fock term since we have only contact 
interactions. Furthermore, the actual expressions are quite difficult to evaluate in a finite 
system. 

As a first approximation, we use the LDA, which may be considered the lowest-order 
term in a derivative expansion of the energy functional. The idea is to expand around the 
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uniform system where (as was shown in Section [II B|) the energy functional can be written 
as an explicit function of p. The LDA prescription is 

^'°^[p(x)] = j rf=^xf2+(po)|po^P(x) , (71) 

where £2+ is the EFT energy density of the uniform system, including terms at second order 
and higher. 

By combining results from Sees. P B| and |T|, we may write the exchange-correlation 
energy, to third order (NNLO) in the EFT expansion, as 

^."°^[P(X)] = j d'^ {^2(P(X)) + ^3(P(X)) + ■ • ■} 

+ (^2 al r, + 63 4 + h al) ^ j rf^x [pi^)f' + ■■■ , (72) 
where the dimensionless bi are 

^.^^(.-1)(-) (ll-21n2), 
1 . /67r2^'/' 



1 , /67r2^'/' 



OTT 

.2\ 5/3 



(0.0755 ((7 - 1) + 0.0574 {g - l)ig - 3)) . 



(73) 



In order to solve for the orbitals in Eq. ( ^9]) and to calculate the energy we need the expression 
for Jo(x). In the LDA, this is simple since 

^"^""^ = 6K^) y'^P^ + E rap] j = (^hfH + . (74) 

To NNLO we find: 

^o(x) = ^ p(x) - I 6, ^ [p(x)]^/3 - I ih, al r, + 63 + K af) ^[p(x)]^/^ • 

(75) 

A convenient expression for the total binding energy (through NNLO) follows by substituting 
for Jo(x) and combining terms: 

^[p(x)] = /^'"[^(")]'-^^^ /^^x[p(x)r 

- ^ {h2 al + 63 + 64 a^) ^ / c^^x [p(x)]«/=' . (76) 

In the numerical calculations given below, the noninteracting case (Co = 0) uses Jo(x) = 
and the first term in Eq. (|76|), LO uses the first term in Eq. (^) and the first two terms in 
Eq. (^, and so on for NLO and NNLO. 
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B. Kohn-Sham Self-Consistent Procedure 



Here we describe the numerical procedure used to find the Kohn-Sham orbitals. We 
assume closed shells, so the density and potentials are functions only of the radial coordinate 
r = |x|. (This restriction is straightforward to relax.) The Kohn-Sham iteration procedure 
is as follows: 

1. Guess an initial density profile p(r). If the system is particularly nonlinear this initial 
choice may be critical. More generally, there may be a metastable state (e.g., if the 
system ultimately collapses). For nuclear systems, the experience is that a crude 
caricature of the true density (such as a Woods-Saxon shape) is adequate. In the 
present case, the non- interacting harmonic oscillator density is sufficient. 

2. Evaluate the local single-particle potential 

Vs[p{r)] = Vs{r) =v{r) - Jo(r) (77) 

using Eqs. (^) at the chosen level of approximation (e.g., NLO). Beyond the LDA Vg 
could be a functional not just of the density but of the Kohn-Sham orbitals individually. 

3. Solve the Schrodinger equation for the lowest A states (including degeneracies), to find 
a set of orbitals and Kohn-Sham eigenvalues {ipa,£a}- 

+ Vsir)] Lfai^) = £a<^„(x) . (78) 



' 2M 

4. Compute a new density from the orbitals: 

A 



a=l 



All other ground-state observables are functionals of 

5. Repeat steps 2.-4. until changes are acceptably small ("self-consistency"). In practice, 
the changes in the density are "damped" by using a weighted average of the densities 
from the {n — l)th and nth iterations: 

p(r)=/3p„_i(r) + (l-/3)p„(r) , (80) 

with < /5 < 1. 

This procedure has been implemented for dilute fermions in a trap using computer codes 
written both in C and in Mathematica. Two methods for carrying out step 3. were tested. 
The Kohn-Sham single-particle equations are solved in one approach by direct integration of 



the differential equations via the Numerov method |4^ and in the other approach by diago- 
nalization of the single-particle Hamiltonian in a truncated basis of unperturbed harmonic 
oscillator wavef unctions. The same results are obtained to high accuracy. For closed shells, 
either method is efficient and easy to code. 
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C. Fermions in a Harmonic Trap 



The interaction through NNLO is specified in terms of the three effective range parameters 
Qg, Vg, and Op. For the numerical calculations presented here, we consider two cases in the 
dilute limit, -C {r^, ttp} ~ with both signs for a^, and hard sphere repulsion with radius 
R, in which case ag = ap = R and Vg = 2R/3. 

Lengths are measured in units of the oscillator parameter b = a/ h/Muj, masses in terms 
of the fermion mass M, and h = 1. In these units, fiw for the oscillator is unity and the 
Fermi energy of a non-interacting gas with filled shells up to Np is Ep = {Np + 3/2). The 
total number of fermions A is related to Np by 



A = ^{Np + l){Np + 2){Np + 3) 
o 



(81) 



Since we have only considered spin-independent interactions, our results are independent of 
whether the spin degeneracy g actually originates from spin, isospin, or some fiavor index. 

With interactions included, single-particle states are labeled by a radial quantum number 
n, an orbital angular momentum / with ^-component mi, and the spin projection. The radial 
functions depend only on n and /, so the degeneracy of each level is 5^ x {21 + 1). Excluding 
spin, the solutions are of the form 



where the radial function Uniir) satisfies 



Uni{r) 



^Imi (^) 



(82) 



Id? , , 1(1 + 1] 



The Unz's are normalized according to 



Unlir) = EnlUnlir) 



\ur,i{r)? dr = 1 



(83) 



(84) 



Thus the density is given by: 



Uiir 



nl 



(85) 



The interactions are sufficiently weak that the occupied states are in one-to-one correspon- 
dence with those occupied in the non-interacting harmonic oscillator potential. 

To check our numerical calculations, we first compared to density distributions at zero 
temperature from Ref. 
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which used a contact interaction with strength corresponding 
to = —0.16 in our units. In that work, non-interacting and mean-field (Hartree-Fock) 
results for the normal state (for possible comparison to superfiuid solutions) were presented 
for systems with 240 and 330 atoms. The corresponding densities in our EFT expansion 
are the curves labeled "Cq = exact" and "Kohn-Sham LO" in Figs. M and O. We also 



show densities for the same systems but with ag = +0.16 in Figs. ^ and |T3|. As one might 
expect, the attractive interaction pulls the density in while the repulsive interaction pushes 
it out relative to the non-interacting density. Values for the energy per particle, the average 
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FIG. 8: Kohn-Sham approximations (see text) for a dilute gas of fcrmions in a harmonic trap with 
degeneracy g = 2 fihed up to Np = 7, which impUes there are 240 particles in the trap. The 
scattering length is Qg = —0.16 and the other effective range parameters are set to zero. 




FIG. 9: Kohn-Sham approximations (see text) for a dilute gas of fermions in a harmonic trap with 
degeneracy g = 2 filled up to Np = 7, which implies there are 240 particles in the trap. The 
scattering length is Qg = +0.16 and the other effective range parameters are set to zero. 
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FIG. 10: Kohn-Sham approximations (see text) for a dilute gas of fermions in a harmonic trap 
with degeneracy g = 4: fiUed up to Np = 4, which imphes there are 140 particles in the trap. The 
scattering length is Qg = —0.10 and the other effective range parameters are set to zero. 




FIG. 11: Kohn-Sham approximations (see text) for a dilute gas of fermions in a harmonic trap 
with degeneracy g = 4 filled up to Np = 4, which implies there are 140 particles in the trap. The 
scattering length is Qg = +0.10 and the other effective range parameters are set to zero. 
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.... Cq = exact 
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FIG. 12: Kohn-Sham approximations (see text) for a dilute gas of fermions in a harmonic trap 
with degeneracy g = 2 fiUed up to Np = 8, which imphes there are 330 particles in the trap. The 
scattering length is = —0.16 and the other effective range parameters are set to zero. 
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FIG. 13: Kohn-Sham approximations (see text) for a dilute gas of fermions in a harmonic trap 
with degeneracy g = 2 filled up to Np = 8, which implies there are 330 particles in the trap. The 
scattering length is Qg = -1-0.16 and the other effective range parameters are set to zero. 
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FIG. 14: Kohn-Sham approximations (see text) for a dilute gas of fermions in a harmonic trap 
with degeneracy g = 2 fiUed up to Np = 7, which imphes there are 240 particles in the trap. The 
scattering lengths are = Op = 0.16 and the effective range is = 2as/3 (hard sphere repulsion). 
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FIG. 15: Kohn-Sham approximations (see text) for a dilute gas of fermions in a harmonic trap 
with degeneracy g = 4 filled up to Np = 4, which implies there are 140 particles in the trap. The 
scattering lengths are Qg = Up = 0.10 and the effective range is rg = 2as/S (hard sphere repulsion). 
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TABLE I: Energies per particle, averages of the local Fermi momentum kp, and rms radii for a 
variety of different parameters and particle numbers for a dilute Fermi gas in a harmonic trap. See 
the text for a description of units. The effective range and p-wave scattering length are set to zero, 
= Up = 0, except for the two lines where has an asterisk, in which case they are given by 
Us = Up = 3rs/2. 
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Np 




as 


E/A 




V(r2) 


approximation 


2 


7 


240 





6.75 


3.27 


2.60 


Co = exact 


2 


7 


240 


-0.16 


5.98 


3.61 


2.35 


KS LO 


2 


7 


240 


-0.16 


6.25 


3.44 


2.47 


KS NLO (LDA) 


2 


7 


240 


-0.16 


6.23 


3.46 


2.46 


KS NNLO (LDA) 


2 


7 


240 


0.16 


7.36 


3.08 


2.76 


KS LO 


2 


7 


240 


0.16 


7.51 


3.03 


2.81 


KS NLO (LDA) 


2 


7 


240 


0.16 


7.52 


3.02 


2.82 


KS NNLO (LDA) 


2 


7 


240 


0.16* 


7.66 


2.97 


2.87 


KS NNLO (LDA) 


4 


4 


140 




4.50 


2.66 


2.12 


Co = exact 


4 


4 


140 


-0.10 


3.62 


3.27 


1.72 


KS LO 


4 


4 


140 


-0.10 


3.83 


3.01 


1.87 


KS NLO (LDA) 


4 


4 


140 


-0.10 


3.75 


3.12 


1.81 


KS NNLO (LDA) 


4 


4 


140 


0.10 


5.09 


2.44 


2.31 


KS LO 


4 


4 


140 


0.10 


5.16 


2.41 


2.34 


KS NLO (LDA) 


4 


4 


140 


0.10 


5.18 


2.40 


2.35 


KS NNLO (LDA) 


4 


4 


140 


0.10* 


5.20 


2.39 


2.36 


KS NNLO (LDA) 



Fermi momentum, and the rms radius for each calculation are given in Table |. Averages 
are defined as 



(/(x)) = l /rf=^x/(x)p(x) 



(86) 



In each of these figures, the other effective range parameters, and a^, are set to zero. 
We have also shown results for a representative system with spin degeneracy g = 4: with 
weaker attractive (a^ = —0.10) and repulsive {as = +0.10) interactions in Figs. |10| and 
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To illustrate the impact of and Up on the NNLO results, we have also included two 
calculations where the underlying interaction is hard-sphere repulsion with R = 0.16 {g = 2) 
and R = 0.10 {g = 4) in Figs. |I| and [T|. 

From Figs. ^ and ^, we see what seems to be good convergence of the density by NNLO. 
Note that the LO results are not well converged in either case. From Table |, we find that 
the average Fermi momentum 



67r^p(x) 



1/3 



P(x) 



(87) 



is such that the expansion parameter in an LDA sense, (/cf)^*; is equal to or greater than 
one-half, which implies poor convergence. In fact, the convergence is misleading, as seen 
by comparison to the hard-sphere repulsive case in Fig. |14|. (See the discussion of energy 
convergence in Sec. |IV F| .) 
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The systems with (? = 4 have {kp)as ~ 1/4-1/3, and the convergence of even the hard- 
sphere case is good. We note that the instabihty for large g discussed in Ref. is mani- 
fested as a collapse of the density iteration by iteration when the stability criterion is violated 
by the interior density. 



D. Comparison to Thomas- Fermi 

In this section, we compare the LDA Kohn-Sham (KS) results to those from a Thomas- 
Fermi (TF) approximation. By Thomas-Fermi, we mean that the entire energy, including 
the kinetic energy, is calculated in a local density approximation. In particular, the Thomas- 
Fermi kinetic energy functional is 

rather than being computed as the sum of the kinetic energy of Kohn-Sham orbitals. We 
define the Thomas-Fermi potential energy functional to have the same form as the Kohn- 
Sham potential energy functional in the LDA. Historically, the Thomas-Fermi approach, as 
applied to atoms and molecules, was the first attempt to use the (electron) density as the 
basic variable rather than solving for the wavefunction. In its original form, only the Hartree 
term and the external nuclear-electron attractive potential were included in the potential 
energy. We will generalize to define LO, NLO, and NNLO Thomas-Fermi approximations. 

The idea behind Thomas-Fermi is that in each volume element dV we have chemical 
equilibrium, with chemical potential /i. Each volume element is labeled by r (we assume 
spherical symmetry for convenience). Local eigenvalues Ek{r) are computed at each r as if 
in a uniform system: 

Ek{r) = ^ + vs{r) , (89) 

and levels are filled until 

^fcp(r)=/i, (90) 

which defines the local Fermi momentum k-F{r). This is turn defines the Thomas- Fermi 
density pTF{r,fi) for a given chemical potential fi as: 

p..(r.^) = |6?P'''l^-''-Ml'"''f''>''-M . (91) 
\ i^^^<v,^r) 

where the potential Vs{r) includes the external trap potential and the LDA Kohn-Sham 
potential [see Eq. (^)]. The conventional TF procedure combines this equation with an 
equation for the potential (e.g., a Poisson equation in the Coulomb case) by substituting for 
the potential. Here we solve it in a two-step process closely analogous to the Kohn-Sham 
solution procedure. 

For a given choice of fi, the number of fermions Atf is given by: 

POO 

Atf (/i) = 47r / p^p {r, p) dr , (92) 
Jo 
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TABLE II: Comparisons of energies per particle, averages of the local Fermi momentum kp, and 
rms radii between Thomas-Fermi and Kohn-Sham treatments of a dilute Fermi gas in a harmonic 
trap. See the text for a description of units. In all cases, = Op = 0. 
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A 


as 


E/A 


(kp) 


V(r2) 


approximation 


2 


2 


20 


— 


3.00 


2.15 


1.73 


Co = exact 


2 


2 


20 


— 


2.94 


2.17 


1.71 


Co = TF 


2 


2 


20 


-0.16 


2.83 


2.24 


1.66 


KS NNLO (LDA) 


2 


2 


20 


-0.16 


2.77 


2.26 


1.64 


TF NNLO (LDA) 


2 


2 


20 


0.16 


3.22 


2.04 


1.83 


KS NNLO (LDA) 


2 


2 


20 


0.16 


3.15 


2.06 


1.81 


TF NNLO (LDA) 


2 


7 


240 




6.75 


3.27 


2.60 


Co = exact 


2 


7 


240 




6.72 


3.29 


2.59 


Co = TF 


2 


7 


240 


-0.16 


6.23 


3.46 


2.46 


KS NNLO (LDA) 


2 


7 


240 


-0.16 


6.20 


3.47 


2.45 


TF NNLO (LDA) 


2 


7 


240 


0.16 


7.52 


3.02 


2.82 


KS NNLO (LDA) 


2 


7 


240 


0.16 


7.49 


3.03 


2.81 


TF NNLO (LDA) 



where we've assumed a spherically symmetric distribution. We find the correct value of fx 
for an input value of A using a root finding program, which finds the zero of 

f{ft)^A-ATF{fi) (93) 
in the interval /imin < < yUmax- The procedure to find the self-consistent ptf is iterative: 

1. Guess an initial pTpir) (for example, the unperturbed density); 

2. given pTpif'), compute Vs{r) for the noninteracting, LO, NLO, or NNLO case; 

3. using this Vs{r) in Eq. (PT|), find /i so that /(/i) = 0; 

4. with the new value of /i, calculate a new pxpi^, fi) from Eq. (^Tj); 

5. return to step 2., continuing until ptf and p do not change within a prescribed toler- 
ance. 

This procedure is rather simple and scales very well with the number of particles. There 
are deficiencies to the Thomas-Fermi approach, however, for many problems of interest. 
These deficiencies are most evident for the Coulomb problem, where it is proved that 



molecules do not bind (the separate atoms always have lower energy) [0. Gradients ex- 
pansions can improve the approximation but have not had a quantitative impact compared 
to Kohn-Sham approaches. 

The kinetic energy contribution is a leading source of nonlocality in the energy functional. 
One of the chief virtues of the Kohn-Sham approach is that it treats this component much 
more effectively than a derivative expansion would. The comparison with Thomas-Fermi 
calculations highlights this difference. The most visible consequences are the absence of 
shell structure in ground state densities in the Thomas-Fermi approach and the incorrect 
treatment of the low-density asymptotic region (dominated by the last Kohn-Sham orbitals). 
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In contrast, the last Kohn-Sham orbital has the correct energy (the ionization energies in 
the exact and Kohn-Sham systems are equal), so the exponential tail of the distribution 
is correct. These deficiencies of Thomas- Fermi are relevant for the calculation of nuclear 
charge and matter densities. 

In Figs. |TB] and Thomas-Fermi and Kohn-Sham densities are compared for a small 
system of 20 atoms (the other parameters are given in the figures). Both non-interacting and 
NNLO curves are shown in each figure. The shell structure is evident in the non-interacting 
density and is only slightly damped by the interactions in the Kohn-Sham approach. In 
contrast, the Thomas- Fermi curves are featureless. The same comparisons but with an 



order of magnitude more atoms {A = 240) are shown in Figs. ^ and |T9|. In these case the 



differences are much smaller. For a small number of particles, the Kohn-Sham procedure is 
a comparable computation to Thomas-Fermi. With thousands of atoms, the Thomas-Fermi 
approximation should be accurate and efficient. 

In Table |l|, the energies per particle and other properties are compared for Kohn-Sham 
and Thomas-Fermi solutions. With A = 20 atoms, the energy is reproduced at about the 2% 
level and the average Fermi momentum at the 1% level. With A = 240 atoms, the energy is 
reproduced better than 1% and the average Fermi momentum to a third of a percent. Thus 
the Thomas-Fermi density will be quite adequate in making error estimates, which we turn 
to next. 



E. Power Counting and Convergence 

A major motivation for the use of effective field theory for many-body problems is the 
promise of error estimates. We face two challenges in making such estimates for a Kohn- 
Sham DFT. There are usually new low-energy constants (LEC's) at each successive order in 
the EFT expansion. We first need to estimate the size of the LEC's in the omitted orders. 
Second, we need to estimate their numerical impact on the functional and subsequently on 
observables in finite systems. In the present case, where we have a perturbative low-density 
expansion, we can do both directly. 

Naive dimensional analysis, or NDA, is frequently used to estimate unknown LEC's in an 
EFT Lagrangian. The idea is to identify the relevant dimensional momentum and mass scales 
in the problem and to rescale each term in the Lagrangian as an appropriate combination 
of scales times a dimensionless coefficient. An estimate of the truncation error follows by 
assuming the dimensionless coefficient is order unity (e.g., from 1/3 to 3). 

In Ref. NDA appropriate to low-energy chiral effective field theories of QCD was 



applied to covariant energy functionals for nuclei. While these functionals were viewed as 
approximate Kohn-Sham energy functionals, their form was directly derived from a chiral 
Lagrangian in the Hartree approximation. That is, there was a one-to-one correspondence 
between terms in the Lagrangian and terms in the corresponding energy functional. As a 
result, estimates of coefficients in the Lagrangian were immediately translated into error 
estimates in the functional, which led to estimates in specific finite nuclei through local 



density approximations |^ 



In the present case, the Hartree (actually Hartree-Fock) terms again provide an imme- 
diate connection between NDA estimates in the Lagrangian and estimates of the energy 
per particle in terms of the fermion density. Using LDA estimates of the average density, 
we get energy estimates for specific systems of trapped atoms. EFT power counting then 
provides estimates for the higher-order terms using the normalization established by the 
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FIG. 16: Thomas-Fermi and Kohn-Sham approximations for dilute gas of fermions in a harmonic 
trap with degeneracy g = 2 filled up to Np = 2, which implies there are 20 particles in the trap. 
The scattering length is Ug = —0.16 and the other effective range parameters are set to zero. 



C„ = exact 




FIG. 17: Thomas-Fermi and Kohn-Sham approximations for dilute gas of fermions in a harmonic 
trap with degeneracy g = 2 filled up to Np = 2, which implies there are 20 particles in the trap. 
The scattering length is Ug = -1-0.16 and the other effective range parameters are set to zero. 
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FIG. 18: Thomas-Fermi and Kohn-Sham approximations for dilute gas of fermions in a harmonic 
trap with degeneracy g = 2 filled up to Np = 7, which imphes there are 240 particles in the trap. 
The scattering length is Ug = —0.16 and the other effective range parameters are set to zero. 




r/b 

FIG. 19: Thomas-Fermi and Kohn-Sham approximations for dilute gas of fermions in a harmonic 
trap with degeneracy g = 2 filled up to Np = 7, which implies there are 240 particles in the trap. 
The scattering length is Ug = -1-0.16 and the other effective range parameters are set to zero. 
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FIG. 20: Contributions to the energy per particle for three example systems of atoms in a harmonic 
trap. The actual contribution in each case is shown as a round symbol. The square symbols are 
estimates based on naive dimensional analysis (see text), with error bars marking a range from 1/2 
to 2 in the dimensionless coefficients. 



Hartree terms. For a purely short-ranged interaction, the NDA is simple: the estimate of 
the Hartree- Fock energy contribution from a given term in the Lagrangian is found by replac- 
ing ip'^tp by the average density (and including an appropriate spin factor). Any reasonable 
density (e.g., the Thomas- Fermi result) can be used to find the average density, since we 
already allow for a much larger uncertainty in the coefficients. Local density estimates for 
derivatives of densities are also sufficiently accurate. 

In Fig. [2^, the contributions to the energy per particle from LO, NLO, and NNLO 
terms in the energy per particle are shown as round symbols for three of the systems from 
Table |. So, for example, the NLO contribution is from |-Enlo — -^loIM- The square symbol 
denote estimates based on naive dimensional analysis, with error bars indicating a 1/2 to 
2 uncertainty in the coefficients. In particular, the LO contribution is an estimate of the 
Hartree-Fock term and the NLO and NNLO estimates were found by multiplying the LO 
result by {kp)as and {{kp)asY, respectively. With one exception, the NDA estimates are 
good predictors of the actual contribution. Therefore, we can use these estimates to predict 
the uncertainty in the energy per particle from higher orders. The exception is the NNLO 
estimate for the g = 2 system, which greatly overestimates the actual contribution at that 
order. The reason is evident from the 64 coefficient in Eq. (^), which determines the size 
of this contribution. Each of the two terms in 64 are the size expected from NDA, but 
for g = 2 they happen to largely cancel. There is always this possibility of unnaturally 
small coefficients (and subsequent contributions) because of accidental cancellations. If we 
compare instead the NNLO energy from the hard-sphere-repulsion case, which has additional 
contributions at NNLO, we find that it is close to the NDA estimate. 
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V. SUMMARY AND CONCLUSIONS 



In this paper, we construct a Kohn-Sham density functional for a confined, dilute Fermi 
gas as an order-by-order effective field theory (EFT) expansion. The starting point is a 
generating functional with a source J{x) coupled to the composite density operator tp'^tp. A 
functional Legendre transformation with respect to the source leads to an effective action of 
the density. In general we might expect complications with such an effective action because 



of the renormalization of composite operators [42, 43 1, but these are avoided in the present 



case because the fermion number is a conserved charge ^ . 

A conventional effective action is a functional of the expectation value of elementary fields 
in the Lagrangian. It defines a classical field theory that contains all of the quantum effects 
of the interacting quantum field theory associated with the original Lagrangian. That is, one 
reproduces results of the full field theory from tree level calculations based on the propagators 
and vertices of the effective action [T^, |T^, |TB|. In the present case, the effective action of 
the density can be used to calculate the ground state energy, including all correlations, with 
what looks like a Hartree calculation (which is the many-body equivalent of tree level). This 
is density functional theory. 



The calculation is carried out by adapting the inversion method proposed in Refs. |22, 23 



and the Kohn-Sham procedure of Ref. to an EFT treatment of the dilute Fermi gas. 



Instead of organizing the perturbative inversion in terms of a coupling constant (e.g., the 
electron charge squared), we use the EFT expansion parameter for the natural dilute sys- 
tem, which is 1/A, where A is the resolution scale. In the uniform system, the ultimate 
dimensionless expansion parameters are products of the Fermi momentum kp and parame- 
ters of the effective range expansion (e.g., the scattering length a^, which is of order 1/A for 
a natural system). In a finite system, the density- weighted average of a local Fermi momen- 
tum times the effective range parameters controls one expansion, with another expansion 
involving gradients of the density. 

In the present work, we used the local density approximation (LDA), which means that 
only the first expansion was tested. The observed convergence of the density and energy in 
sample systems confirms this expansion for a finite system. An error plot of contributions to 
the energy per particle versus the order of the calculation shows that we can reliably estimate 
the truncation error in a finite system. The next step is to develop a derivative expansion 
and test it for convergence. Work is in progress to adapt to the present problem the methods 
that have been used to derive derivative expansions for one-loop effective actions. 

The EFT expansion for the uniform system using dimensional regularization and minimal 
subtraction (DR/MS) offered many simplifications over conventional treatments |34|. The 
derivative expansion approach would preserve all of these advantages. In particular, renor- 
malization of the derivative expansion coefficient functions would be carried out entirely 
in the uniform system using DR/MS, so that a given diagram contributes to only a single 
prescribed order in the expansion. 

An alternative to the inversion method that might be preferred for some systems is 
to introduce an auxiliary field (or fields) and to carry out the Legendre transformation 



conventionally with respect to that field p2l 23]. This could be applied to the large-N EFT 



expansion discussed in Ref. [46|. The construction of Kohn-Sham DFT in the auxiliary field 



formulation is discussed in Refs. \U7[ BS 



Many of the phenomena of greatest interest in experiments on trapped fermion atoms 
involve tuning the system so that the s-wave scattering length is large. For example, one 



34 



would like to study superfiuid transitions at low temperature W7[ 091 pO|. systematic solu- 



tion to the large scattering length problem at finite density is not yet available, even for a 
uniform system. Furthermore, with any attractive interaction we expect a pairing instabil- 
ity. Extending the DFT/EFT expansion procedure to include pairing will be an important 
challenge. The inversion method has been applied to conventional BCS superconductivity 
in Ref. |^ using a source coupled to the pair creation and destruction operator. Work to 
adapt this procedure to the EFT is in progress. 

These extensions will be directly relevant for nuclear applications as well. In addition, to 
adapt the density functional procedure to chiral effective field theories with explicit pions, 
we will need to extend the discussion to include long-range forces ("long range" means 
compared to 1/A). Other extensions include DFT with alternative source terms (e.g., spin- 
density DFT) and time-dependent DFT in the EFT formalism. Studies of each of these 
extensions are in progress. 
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